# loading packages/libraries
library(sf)
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(here)
## here() starts at D:/CASA/GIS/gis_code/Prac5
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tmap)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(tmaptools)
library(grid)
library(leafpop)
library(leaflet)
# reading in the London Boroughs data
Londonborough <- st_read(here("prac5_data",
"statistical-gis-boundaries-london",
"ESRI",
"London_Borough_Excluding_MHW.shp")) %>%
select(c(-SUB_2009,-SUB_2006)) %>%
st_transform(., 27700) # transforming the CRS of this data set
## Reading layer `London_Borough_Excluding_MHW' from data source
## `D:\CASA\GIS\gis_code\Prac5\prac5_data\statistical-gis-boundaries-london\ESRI\London_Borough_Excluding_MHW.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
## Projected CRS: OSGB36 / British National Grid
# reading the OSM data
OSM <- st_read(here("prac5_data",
"greater-london-latest-free.shp",
"gis_osm_pois_a_free_1.shp")) %>%
st_transform(., 27700) %>%
#select hotels only
dplyr::filter(fclass == 'hotel')
## Reading layer `gis_osm_pois_a_free_1' from data source
## `D:\CASA\GIS\gis_code\Prac5\prac5_data\greater-london-latest-free.shp\gis_osm_pois_a_free_1.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 48312 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -0.5108706 ymin: 51.28109 xmax: 0.322123 ymax: 51.6926
## Geodetic CRS: WGS 84
# a spatial join example for these two datasets
join_example <- st_join(Londonborough, OSM)
head(join_example)
## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 516362.6 ymin: 159907.4 xmax: 522654 ymax: 172367
## Projected CRS: OSGB36 / British National Grid
## NAME GSS_CODE HECTARES NONLD_AREA ONS_INNER osm_id code
## 1 Kingston upon Thames E09000021 3726.117 0 F 28055159 2401
## 1.1 Kingston upon Thames E09000021 3726.117 0 F 161631526 2401
## 1.2 Kingston upon Thames E09000021 3726.117 0 F 421483573 2401
## 1.3 Kingston upon Thames E09000021 3726.117 0 F 421506999 2401
## 1.4 Kingston upon Thames E09000021 3726.117 0 F 502772206 2401
## 1.5 Kingston upon Thames E09000021 3726.117 0 F 893668667 2401
## fclass name
## 1 hotel Premier Inn
## 1.1 hotel Chessington Safari Hotel
## 1.2 hotel Premier Inn Chessington
## 1.3 hotel Chessington Azteca Hotel
## 1.4 hotel Warren House
## 1.5 hotel DoubleTree by Hilton Kingston Upon Thames
## geometry
## 1 MULTIPOLYGON (((516401.6 16...
## 1.1 MULTIPOLYGON (((516401.6 16...
## 1.2 MULTIPOLYGON (((516401.6 16...
## 1.3 MULTIPOLYGON (((516401.6 16...
## 1.4 MULTIPOLYGON (((516401.6 16...
## 1.5 MULTIPOLYGON (((516401.6 16...
# won't run because disabled by eval=FALSE, and because I have done this above only here for notes purposes as we shift into Static Map
#Load all our data
library(sf)
library(tmap)
library(tmaptools)
library(tidyverse)
library(here)
# read in all the spatial data and
# reproject it
OSM <- st_read(here::here("prac5_data",
"greater-london-latest-free.shp",
"gis_osm_pois_a_free_1.shp")) %>%
st_transform(., 27700) %>%
#select hotels only
filter(fclass == 'hotel')
##London Borough data is already in 277000
Londonborough <- st_read(here::here("Prac1_data",
"statistical-gis-boundaries-london",
"ESRI",
"London_Borough_Excluding_MHW.shp"))%>%
st_transform(., 27700)
# loading more data
Worldcities <- st_read(here::here("prac5_data",
"World_Cities",
"World_Cities.shp")) %>%
st_transform(., 27700)
## Reading layer `World_Cities' from data source
## `D:\CASA\GIS\gis_code\Prac5\prac5_data\World_Cities\World_Cities.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2540 features and 14 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -19609100 ymin: -7321602 xmax: 19950890 ymax: 14476660
## Projected CRS: WGS 84 / Pseudo-Mercator
UK_outline <- st_read(here::here("prac5_data",
"gadm41_GBR_shp",
"gadm41_GBR_0.shp")) %>%
st_transform(., 27700)
## Reading layer `gadm41_GBR_0' from data source
## `D:\CASA\GIS\gis_code\Prac5\prac5_data\gadm41_GBR_shp\gadm41_GBR_0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -8.649996 ymin: 49.86531 xmax: 1.764393 ymax: 60.84548
## Geodetic CRS: WGS 84
Example code to unzip a .gz file
# read in the .csv
# and make it into spatial data
Airbnb <- read_csv(here("prac5_data",
"listings.csv")) %>%
# longitude is considered x value here, latitude is y
st_as_sf(.,
coords = c("longitude", "latitude"),
crs = 4326) %>%
st_transform(., 27700)%>%
#select entire places that are available all year
filter(room_type == 'Entire home/apt' & availability_365 =='365')
## Rows: 96182 Columns: 18
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): name, host_name, neighbourhood, room_type
## dbl (11): id, host_id, latitude, longitude, price, minimum_nights, number_o...
## lgl (2): neighbourhood_group, license
## date (1): last_review
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Lets make a function to join our data. It will aggregate counts of the number of times a borough appears
# make a function for the join
# functions are covered in practical 7
# but see if you can work out what is going on
# hint all you have to do is replace data1 and data2
# with the data you want to use
Joinfun <- function(data1, data2){
output<- data1%>%
st_join(data2,.) %>% # notice the second data set is on the left
add_count(GSS_CODE, name="hotels_in_borough") # function creating counts of GSS_CODE and storing in column "hotels_in_borough"
return(output)
}
Use the new function to join our data sets
# use the function for hotels
Hotels <- Joinfun(OSM, Londonborough)
# then for airbnb
# this is incorrect - change to airbnb2 to look at result
Airbnb <- Joinfun(Airbnb, Londonborough)
Worldcities2 <- Worldcities %>%
filter(CNTRY_NAME=='United Kingdom'&
Worldcities$CITY_NAME=='Birmingham'|
Worldcities$CITY_NAME=='London'|
Worldcities$CITY_NAME=='Edinburgh')
newbb <- c(xmin=-296000, ymin=5408, xmax=655696, ymax=1000000) # setting a bounding box
UK_outlinecrop <- UK_outline$geometry %>%
st_crop(., newbb) # doing a spatial crop using the new bounding box
Hotels <- Hotels %>%
#at the moment each hotel is a row for the borough
#we just want one row that has number of airbnbs
group_by(., GSS_CODE, NAME) %>%
summarise(`Accomodation count` = unique(hotels_in_borough))
## `summarise()` has grouped output by 'GSS_CODE'. You can override using the
## `.groups` argument.
Airbnb <- Airbnb %>%
group_by(., GSS_CODE, NAME)%>%
summarise(`Accomodation count` = unique(hotels_in_borough))
## `summarise()` has grouped output by 'GSS_CODE'. You can override using the
## `.groups` argument.
Hotels %>%
filter(NAME=="Sutton")
## Simple feature collection with 1 feature and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 522223.3 ymin: 159658.1 xmax: 531245.8 ymax: 167622.5
## Projected CRS: OSGB36 / British National Grid
## # A tibble: 1 × 4
## # Groups: GSS_CODE [1]
## GSS_CODE NAME `Accomodation count` geometry
## * <chr> <chr> <int> <MULTIPOLYGON [m]>
## 1 E09000029 Sutton 1 (((528552.3 159658.1, 528399.7 159928.8…
tmap_mode("plot")
## tmap mode set to plotting
# set the breaks for our mapped data
breaks = c(0, 5, 12, 26, 57, 286)
# plot each map
tm1 <- tm_shape(Hotels) +
tm_polygons("Accomodation count",
breaks=breaks,
palette="PuBu")+
tm_legend(show=FALSE)+
tm_layout(frame=FALSE)+
tm_credits("(a)", position=c(0,0.85), size=1.5)
tm2 <- tm_shape(Airbnb) +
tm_polygons("Accomodation count",
breaks=breaks,
palette="PuBu") +
tm_legend(show=FALSE)+
tm_layout(frame=FALSE)+
tm_credits("(b)", position=c(0,0.85), size=1.5)
tm3 <- tm_shape(UK_outlinecrop)+
tm_polygons(col="darkslategray1")+
tm_layout(frame=FALSE)+
tm_shape(Worldcities2) +
tm_symbols(col = "red", scale = .5)+
tm_text("CITY_NAME", xmod=-1, ymod=-0.5)
legend <- tm_shape(Hotels) +
tm_polygons("Accomodation count",
breaks=breaks,
palette="PuBu") +
tm_scale_bar(position=c(0.2,0.04), text.size=0.6)+
tm_compass(north=0, position=c(0.65,0.6))+
tm_layout(legend.only = TRUE, legend.position=c(0.2,0.25),asp=0.1)+
tm_credits("(c) OpenStreetMap contrbutors and Air b n b", position=c(0.0,0.0))
t=tmap_arrange(tm1, tm2, tm3, legend, ncol=2)
t
## Warning: Values have found that are higher than the highest break
## Warning: Values have found that are higher than the highest break
#library(grid)
# erases the current device or moves to a new page
# probably not needed but makes sure you are plotting on a new page.
grid.newpage()
pushViewport(viewport(layout=grid.layout(2,2)))
print(tm1, vp=viewport(layout.pos.col=1, layout.pos.row=1, height=5))
## Warning: Values have found that are higher than the highest break
print(tm2, vp=viewport(layout.pos.col=2, layout.pos.row=1, height=5))
print(tm3, vp=viewport(layout.pos.col=1, layout.pos.row=2, height=5))
print(legend, vp=viewport(layout.pos.col=2, layout.pos.row=2, height=5))
## Warning: Values have found that are higher than the highest break
# calculating the bounding box for the London Airbnb dataset
Londonbb <- st_bbox(Airbnb,
crs = st_crs(Airbnb))%>%
#we need this to convert it into a class of sf
# otherwise if our bb won't have a class it will just be x and y coordinates for the box
# this makes it into a polygon
st_as_sfc()
main <- tm_shape(Airbnb, bbbox = Londonbb) +
tm_polygons("Accomodation count",
breaks=breaks,
palette="PuBu")+
tm_scale_bar(position = c("left", "bottom"), text.size = .75)+
tm_layout(legend.position = c("right","top"),
legend.text.size=.75,
legend.title.size = 1.1,
frame=FALSE)+
tm_credits("(c) OpenStreetMap contrbutors and Air b n b", position=c(0.0,0.0))+
#tm_text(text = "NAME", size = .5, along.lines =T, remove.overlap=T, auto.placement=F)+
tm_compass(type = "8star", position = c(0.06, 0.1)) +
#bottom left top right <- all round margin
tm_layout(inner.margin=c(0.02,0.02,0.02,0.2))
inset = tm_shape(UK_outlinecrop) +
tm_polygons() +
tm_shape(Londonbb)+
tm_borders(col = "grey40", lwd = 3)+ ##lwd==line width
tm_layout(frame=FALSE,
bg.color = "transparent")+
tm_shape(Worldcities2) +
tm_symbols(col = "red", scale = .5)+
tm_text("CITY_NAME", xmod=-1.5, ymod=-0.5) # positioning of the text with xmod and ymod, -ve means left and down respectively
#library(grid)
main
print(inset, vp = viewport(0.86, 0.29, width = 0.5, height = 0.55))
tmap_save(t, 'hotelsandairbnbR.png')
## Warning: Values have found that are higher than the highest break
## Warning: Values have found that are higher than the highest break
## Map saved to D:\CASA\GIS\gis_code\Prac5\hotelsandairbnbR.png
## Resolution: 2100 by 2100 pixels
## Size: 7 by 7 inches (300 dpi)
#library(grid)
tmap_save(tm = main,insets_tm = inset, insets_vp = viewport(0.86, 0.29, width=.5, height=.55), filename="test.pdf", dpi=600)
## Map saved to D:\CASA\GIS\gis_code\Prac5\test.pdf
## Size: 8.833333 by 5.541667 inches
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(Airbnb) +
tm_polygons("Accomodation count", breaks=breaks)
# library for pop up boxes
#library(leafpop)
#library(leaflet)
#join data
Joined <- Airbnb%>%
st_join(., Hotels, join = st_equals)%>%
dplyr::select(GSS_CODE.x, NAME.x, `Accomodation count.x`, `Accomodation count.y`)%>%
dplyr::rename(`GSS code` =`GSS_CODE.x`,
`Borough` = `NAME.x`,
`Airbnb count` = `Accomodation count.x`,
`Hotel count`= `Accomodation count.y`)%>%
st_transform(., 4326)
#remove the geometry for our pop up boxes to avoid
popupairbnb <-Joined %>%
st_drop_geometry()%>%
dplyr::select(`Airbnb count`, Borough)%>%
popupTable()
## Adding missing grouping variables: `GSS code`
popuphotel <-Joined %>%
st_drop_geometry()%>%
dplyr::select(`Hotel count`, Borough)%>%
popupTable()
## Adding missing grouping variables: `GSS code`
tmap_mode("view")
## tmap mode set to interactive viewing
# set the colour palettes using our previously defined breaks
pal1 <- Joined %>%
colorBin(palette = "YlOrRd", domain=.$`Airbnb count`, bins=breaks)
pal1 <-colorBin(palette = "YlOrRd", domain=Joined$`Airbnb count`, bins=breaks)
pal2 <- Joined %>%
colorBin(palette = "YlOrRd", domain=.$`Hotel count`, bins=breaks)
map<- leaflet(Joined) %>%
#add our polygons, linking to the tables we just made
addPolygons(color="white",
weight = 2,
opacity = 1,
dashArray = "3",
popup = popupairbnb,
fillOpacity = 0.7,
fillColor = ~pal2(`Airbnb count`),
group = "Airbnb")%>%
addPolygons(fillColor = ~pal2(`Hotel count`),
weight = 2,
opacity = 1,
color = "white",
dashArray = "3",
popup = popupairbnb,
fillOpacity = 0.7,group = "Hotels")%>%
#add basemaps
addTiles(group = "OSM (default)") %>%
addProviderTiles(providers$Stadia.StamenToner, group = "Toner") %>%
addProviderTiles(providers$Stadia.StamenTonerLite, group = "Toner Lite") %>%
addProviderTiles(providers$CartoDB.Positron, group = "CartoDB")%>%
# add a legend
addLegend(pal = pal2, values = ~`Hotel count`, group = c("Airbnb","Hotel"),
position ="bottomleft", title = "Accomodation count") %>%
# specify layers control
addLayersControl(
baseGroups = c("OSM (default)", "Toner", "Toner Lite", "CartoDB"),
overlayGroups = c("Airbnb", "Hotels"),
options = layersControlOptions(collapsed = FALSE)
)
## Warning in pal2(`Hotel count`): Some values were outside the color scale and
## will be treated as NA
# plot the map
map
# this is what was explained in Prac6 section: 6.5.4
# attempting a spatial join, by default, it uses st_intersects which is problematic. See below
all_accomodation <- st_join(Hotels,Airbnb)
head(all_accomodation)
## Simple feature collection with 6 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 530967.7 ymin: 180510.3 xmax: 533839.6 ymax: 182206.1
## Projected CRS: OSGB36 / British National Grid
## # A tibble: 6 × 7
## # Groups: GSS_CODE.x [1]
## GSS_CODE.x NAME.x `Accomodation count.x` geometry GSS_CODE.y
## <chr> <chr> <int> <MULTIPOLYGON [m]> <chr>
## 1 E09000001 City o… 23 (((531145.1 180782.1, 53… E09000001
## 2 E09000001 City o… 23 (((531145.1 180782.1, 53… E09000007
## 3 E09000001 City o… 23 (((531145.1 180782.1, 53… E09000012
## 4 E09000001 City o… 23 (((531145.1 180782.1, 53… E09000019
## 5 E09000001 City o… 23 (((531145.1 180782.1, 53… E09000030
## 6 E09000001 City o… 23 (((531145.1 180782.1, 53… E09000033
## # ℹ 2 more variables: NAME.y <chr>, `Accomodation count.y` <int>
# lets correct this by changing the argument to st_equals
all_accomodation2 <- st_join(Hotels, Airbnb, join = st_equals)
head(all_accomodation2)
## Simple feature collection with 6 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 515484.9 ymin: 156480.9 xmax: 554089.2 ymax: 198355.2
## Projected CRS: OSGB36 / British National Grid
## # A tibble: 6 × 7
## # Groups: GSS_CODE.x [6]
## GSS_CODE.x NAME.x `Accomodation count.x` geometry GSS_CODE.y
## <chr> <chr> <int> <MULTIPOLYGON [m]> <chr>
## 1 E09000001 City o… 23 (((531145.1 180782.1, 53… E09000001
## 2 E09000002 Barkin… 9 (((543905.4 183199.1, 54… E09000002
## 3 E09000003 Barnet 14 (((524579.9 198355.2, 52… E09000003
## 4 E09000004 Bexley 6 (((547226.2 181299.3, 54… E09000004
## 5 E09000005 Brent 13 (((525201 182512.6, 5251… E09000005
## 6 E09000006 Bromley 4 (((540373.6 157530.4, 54… E09000006
## # ℹ 2 more variables: NAME.y <chr>, `Accomodation count.y` <int>